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ABSTRACT 


Formulas  from  Lambert’s  problem  in  celestial  mechanics  are  pre¬ 
sented  for  use  in  a  computer  program  to  calculate  velocity  requirements 
for  rapid  intercept  of  an  earth  satellite  by  a  rocket  fired  from  the  earth  or 
from  an  orbit  about  the  earth.  The  rocket  is  assumed  to  receive  an  impul¬ 
sive  velocity  at  launch  and  the  effect  of  impulsive  midcourse  corrections 
are  considered.  The  effects  of  the  earth  rotation  are  included  if  the  rocket 
is  fired  from  the  earth,  but  those  due  to  the  earth’s  atmosphere  are  ignored. 
The  rocket  and  target  satellite  are  assumed  to  be  moving  in  conic  section 
in  the  earth's  gravitational  field. 


Accepted  for  the  Air  Force 

Joseph  R.  Waterman,  Lt.  Col.,  USAF 

Chief,  Lincoln  Laboratory  Project  Office 


TABLE  OF  CONTENTS 


Abstract  iii 

I.  Introduction  1 

II.  Rotation  of  the  Earth  3 

III.  Position  and  Velocity  in  a  Conic  Section  Orbit  given  the 

Orbital  Elements  6 

A.  Parabolic  Motion  13 

B.  Elliptic  Motion  15 

C.  Hyperbolic  Motion  20 

D.  Straight  Line  Motion  24 

IV.  Determination  of  Conic  Section  Orbital  Elements  given  the 

Position  and  Velocity  26 

A.  Parabolic  Motion  27 

B.  Elliptic  Motion  28 

C.  Hyperbolic  Motion  30 

V.  Lambert  Problem  for  the  Time  of  Travel  Between  Two 

Points  in  a  Conic  Section  Orbit  30 

A.  Parabolic  Motion  33 

B.  Elliptic  Motion  37 

C.  Hyperbolic  Motion  42 

VI.  Logical  Flow  of  the  Computer  Program  to  Calculate  Time 

to  Intercept  Versus  Initial  Velocity  47 

References  49 


IV 


VELOCITY  REQUIREMENTS  FOR  RAPID  INTERCEPT 
WITH  MIDCOURSE  CORRECTIONS 


I.  INTRODUCTION 

This  technical  note  documents  formulas  to  be  used  in  a  Fortran 
computer  program  to  calculate  velocity  requirements  for  rapid  intercepts 
of  an  earth  satellite  by  a  rocket  fired  from  the  earth  or  from  an  orbit  about 
the  earth.  The  rocket  is  imagined  to  receive  an  impulsive  velocity  at 
launch.  The  effect  of  impulsive  midcourse  corrections  is  considered.  The 
effects  of  the  earth  rotation  are  included  if  the  rocket  is  fired  from  the 
earth,  but  those  due  to  the  earth's  atmosphere  are  ignored.  The  rocket 
and  the  satellite  are  taken  to  be  moving  in  conic  sections  in  the  —  fj/r 
gravitational  potential  of  the  earth  with  higher  harmonics  ignored,  where  fj. 
is  the  gravitational  constant  times  the  mass  of  the  earth  and  where  r  is 
the  distance  from  the  center  of  the  earth.  Units  used  in  the  program  are 
kilometers  and  kilometers  per  second. 

Input  to  the  program  (with  names  in  the  program  in  capital  letters) 

ar  e : 

(a)  Calendar  date:  IMONTH,  IDAY,  IYEAR 

(b)  Universal  time  of  rocket  firing:  IHR,  IMIN,  SEC 

(c)  Coordinates  of  the  launching  site  (if  the  rocket  is 
launched  from  the  earth) 

p  =  RADIUS  =  distance  from  center  of  earth  (km) 

0  =  LONG  =  longitude  west  of  Greenwich  (deg) 

0  =  LAT  =  north  latitude  (deg)  (1 ) 

(d)  Elliptic  orbital  elements  of  the  target  satellite  in  the 
coordinate  system  referred  to  true  equinox  and  equator 
of  date: 


a  =  A  =  semi-major  axis  (km) 


1 


e 


E 


eccentricity 


I  =  INC  =  inclination  (deg) 

Q  =  ASC  =  right  ascension  of  ascending  node  (deg) 

ou  =  PER  =  argument  of  perigee  (deg) 

M  =  ANOM  =  mean  anomaly  at  time  t  of 

o  y  o 

rocket  firing  (2) 

[For  the  intercept  rocket  orbit,  which  could  be  hyperbolic  or  parabolic  as 
well  as  elliptic,  we  use  internally  in  the  program  instead  of  a  ,  the 

following  elements 


p  =  semi-latus  rectum 

t  =  time  of  perigee  crossing,  ]  (3) 

P 

(e)  Similar  elliptic  orbital  elements  for  inital  parking  orbit 
of  the  rocket  if  it  is  launched  from  orbit, 

(f)  Increments  DVLNCH  of  magnitude  of  velocity  at  launch 
from  an  initial  possible  magnitude  VLNCHO  to  a  final 
possible  magnitude  VLNCHl. 

(g)  Errors  in  the  orbital  elements  of  the  target  satellite 
(presumed  discovered  after  launch)  if  midcourse  cor¬ 
rections  are  considered: 


Aa  =  DA 
Ae  =  DE 
AI  =  DINC 
Afi  =  DASC 


2 


Auu  =  DPER 


AM  =  DANOM 
o 

(h)  Epochs  after  launch  at  which  midcourse  impulsive  firings 
are  possible* 

(i)  Increments  of  magnitude  of  midcourse  velocity  change 
from  an  initial  possible  magnitude  of  change  to  a  final 
possible  magnitude  of  change. 

Output  from  the  program  are: 

(a)  Minimum  time  to  intercept  versus  magnitude  of  impulsive 
velocity  at  launch  at  the  given  increments  of  magnitude  of 
velocity  at  launch; 

(b)  For  each  possible  epoch  of  midcourse  correction  and  for 
each  possible  intercept  orbit  from  launch  to  the  original 
target  orbit,  the  minimum  time  to  intercept  the  new 
target  orbit  versus  magnitude  of  impulsive  velocity  at 
midcourse  correction  at  the  given  increments  of 
magnitude  of  midcourse  velocity  change. 

There  will  be  SC4060  graphical  output  so  that  by  varying  the  input 
one  can  see  the  variation  in  velocity  requirements  versus  time  of  intercept 
as  a  function  of  target  satellite  orbit,  launching  site  or  parking  orbit,  and 
time  of  launch. 

In  the  following  we  rigorously  derive  many  standard  formulas  in  the 
philosophy  that  too  much  documentation  of  a  computer  program  is  better 
than  too  little. 

II.  ROTATION  OF  THE  EARTH 
12  3 

Let  (x  ,  x  ,  x  )  be  a  coordinate  system  referred  to  the  true  equinox 

3 

and  equator  of  date  with  origin  at  the  center  of  the  earth.  The  x  axis 


3 


points  to  the  north  along  the  axis  of  rotation  of  the  earth,  the  x  axis  lies 

2 

in  the  equator  and  points  towards  the  first  point  of  Aries,  and  the  x  axis 
completes  the  right  hand  system. 

Let  (p,  0,  0)  be  the  coordinates  of  the  launch  site  as  defined  in  (1). 
Let  s  be  the  true  sidereal  time,  i.  e.  ,  the  Greenwich  hour  angle  of  the 
first  point  of  Aries.  Then  the  cartesian  coordinates  of  the  launch  site  are 


x^  =  p  cos  0  cos  (s  —  0) 

2 

=  o  cos  0sin  (s  —  9) 


p  sin  0 


(4) 


The  velocity  of  the  rocket  as  it  sits  at  the  launch  site  due  to  the  rotation  of 
the  earth  is 


•  1  2  ds 

x  =  —  x  — 


x  = 

o 

x  dt 

.  2 

1 

ds 

x  = 

o 

X 

dt 

•  3 

x  = 

0 

(5) 

— ♦  — »  12  3 

Let  (e^,  e^,  e^)  be  unit  vectors  in  the  (x  ,  x  ,  x  )  coordinate 

directions.  The  vector  from  the  center  of  the  earth  to  the  launch  site  is 


r  =  x  e,  +  x  e0  +  x  e0  (6) 

o  o  1  o  ^  o  3 


and  the  velocity  of  the  rocket  at  the  launch  site  due  to  the  rotation  of  the 
earth  is 


v 

o 


x 


+  x 


(V) 


4 


How  is  the  sidereal  time  calculated  for  the  given  calendar  date  and 
universal  time  UT  of  launch?  First  we  determine  the  Julian  Day  Number 
JD  from  the  given  month,  day  and  year  of  launch  as  given  in  Ref.  1, 

Table  I,  p.  445.  A  computer  subroutine  can  easily  be  written  to  do  this 

using  the  fact  that  there  are  365  days  in  a  year  except  366  days  in  years 

divisible  by  4  .  The  Julian  Date  at  midnight  beginning  of  day  is  then 

JD  —  0.5  .  The  Greenwich  mean  sidereal  time  s  at  0  universal  time 

o 

on  the  day  of  interest  is 


s  =  6h38m45s.  836  +  8,  640,  1 84®  542T  +  0s.  0929T‘ 
o 


(8) 


where  T  denotes  the  number  of  Julian  centuries  of  36525  days  which,  at 
midnight  beginning  of  day,  have  elapsed  since  mean  noon  on  1900  January  0 
at  the  Greenwich  meridian  (Julian  Date  2415020.  0)  ;  see  Ref.  1  ,  p.  474. 

The  Greenwich  true  sidereal  time  at  a  given  instant  UT  of  universal  time 
on  that  day  is  then 


+  X  UT  +  A*  COS  6 

dt 


(9) 


(1.002737909265  +  0.589  x  10_1°T) 

sidereal  time  seconds  per  universal  time  second  (10) 

and  where  Aijf  is  the  nutation  in  longitude  and  e  the  obliquity  of  the  ecliptic; 
see  Ref.  2,  pp.  75-76.  We  shall  ignore  the  A^  cos  e  term,  whose  largest 
magnitude  is  if 3  .  This  has  the  effect  of  changing  the  launch  site  longitude 
by  this  amount.  For  use  in  (5)  ,  formula  (10)  for  ds/dt  must  be  multi¬ 
plied  by  2tt/86,400  to  convert  to  radians  per  second. 


where 


ds 

dt 
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III. 


POSITION  AND  VELOCITY  IN  A  CONIC  SECTION 
ORBIT  GIVEN  THE  ORBITAL  ELEMENTS 


The  equations  of  motion  of  a  body  of  negligably  small  mass  in  the 
gravitational  field  of  a  spherically  symmetric  earth  are 


d2r 

dt2 


_U£ 

3 

r 


(11) 


where  [i  >  0  is  the  gravitational  constant  times  the  mass  of  the  earth,  r 
is  the  position  vector  of  the  body  and  r  =  |r|  =  (r  •  r)  .  Taking  the 

dot  product  of  both  sides  of  (11)  with  dr/dt  we  obtain 

d^r*  dr  ^  «  dr 

•  a r  "  i  r  •  dt 

dt  r 


which  implies 

1  _d_  /dr  dr \  _  _  1_  ji_  d(r  »  r)  _  ji_  /l\ 
T  dt  \dF  dt)  -  2  r3  dt  dt  \r) 


Thus  we  obtain  the  conservation  of  energy  result 

V2  =  zfe  +  c)  (12) 

where  v  =  dr/dt  and  C  is  a  constant  of  integration.  Taking  the  cross 
product  of  (11)  with  r  we  obtain 

d2r  3 
r  x  — y  -  0 

dt 

which  implies 


6 


Thus  we  obtain  the  conservation  of  angular  momentum  result 


dr  =r 

x  dF  =  c 


(13) 


where  G  is  a  constant  vector  of  integration.  Taking  the  dot  product  of  (13) 
with  r  gives 

G  •  r  =  0  (14) 


which  implies  that  the  motion  is  in  a  plane  perpendicular  to  G  and  that 
there  are  only  two  independent  components  of  G  ,  not  three. 

For  the  moment  we  exclude  the  straight  line  case  and  assume 
G  ^  0  .  Referring  to  Fig.  1,  the  values  of  the  ascending  node  Cl  and  inclin¬ 
ation  I  of  the  orbital  plane  normal  to  G  are  given  by 


G  •  en 

(15) 


(16) 


P 

— »  — »  1  z 

where  G  =  I G  |  and  where  G  is  the  projector  of  G  onto  the  (x  ,  x  ) 

P 

plane, 


cos  I  = 


cos 


Q  - - E- 


G 

G  •  e- 


sin 


G  •  e1 

o  =  — L 


0  £  I  £  180 


0  <;  Cl  <  360 


Gp  =  G  -  (G  •  e3)  e3  (17) 

1  2 

Let  (u  ,  u  )  be  the  coordinates  of  the  body  in  a  coordinate  system 

1  1  2 
with  u  axis  being  the  intersection  of  the  orbit  plane  with  the  (x  ,  x  ) 

2 

plane  and  with  u  axis  in  the  orbit  plane  completing  the  right  hand  system. 
We  define  circular  coordinates  (r,  <jr  )  in  the  orbit  plane  by 


1  ~ 

u  =  r  cos  i)f 

2  ~ 

u  =  r  sin  f  (18) 


7 


Fig.  1.  Relation  between  inclination  I  and  ascending  node  fi 
of  a  plane  and  the  normal  G  to  the  plane. 


8 


Then  (13)  implies 


and  (12)  implies 
2 


+  r 


2  |d  y 
dt 


2  1^-  +  C 


Since 


dr  _  dr_  dj[  _  d_r  _G_  _ _ d_  /jG\ 

dt  ,7  dt  ,7  2  \r/ 

d  ijf  d  i|/  r  d  i|/ 


we  have 


2  ♦  ©* 


2|t  +  c 


£  ©  +  4  +  c) 


Let 


=  +  s. 

G  r 


Q2  =  2  C  +  ^-2 
G 


W e  then  have 


=v 


Q2  -  q2 


which  implies 


(19) 


(20) 


±  (T  -  <n) 


arc  cos 


9 


where  uu  is  a  constant  of  integration.  Since  cos  (+  8) 


cos  (—8)  we  have 


r  =  5  *  h  V2C  +  z 


COS  (lj/  —  U)) 


We  can  have  the  same  curve  with  a  —  sign  as  with  a  +  sign  by  changing 
the  constant  of  integration  uu  by  T7  .  Therefore, 


2-  =  1  4-  e  cos 

r 


(21) 


where 


P  = 


G_ 

4 


>  0  (semi-latus  rectum) 


(22) 


>[- 


2c  ^  c 

e  =  ^ ^ -  +  1  £  0  (eccentricity) 


uu  =  argument  of  perigee 


(23) 

(24) 


—  uu  =  (true  anomaly) 


(25) 


We  have  that  f  =  0  at  perigee,  the  point  of  closest  approach  to  the  earth. 

The  curve  represented  by  (21)  is  one  of  the  following  types: 

0  ^  e  <  1  ellipse 

e  =  1  parabola 

e  >  1  hyperbola  (26) 

with  a  focus  at  the  origin. 

We  define  the  semi-major  axis  a  by 


1  -  e 


(27) 


10 


where 


a  >  0  ellipse 

a  =  0°  parabola 

a  <  0  hyperbola 


(28) 


From  equations  (12),  (22),  (23)  and  (27)  we  have 


2 

v 


(29) 


Let  (y  ,  y  ,  y  )  be  a  coordinate  system  with  origin  at  the  center 

1  3 

of  the  earth  with  y  axis  pointing  towards  perigee,  with  y  axis  normal 

2 

to  the  orbital  plane  along  the  angular  momentum  vector  G  ,  and  with  y 
axis  completing  the  right  hand  system  (see  Fig.  2).  Then  we  have 


x 


j 


3 


E 

k=i 


1,2,3 


k 

y 


3 

E 


j 


i=i 


b.,  xJ 
Jk 


k  =  1,2,3 


(30) 


where  the  orthogonal 


b 

b 

b 

b 

b 


1  1  “ 

12 

13 

21  = 
22  = 


matrix  (b.,  ) 
Jk 

COS  n  COS  CJD  - 

—  cos  0  sin  u) 

sin  n  sin  I 

sinficoscr  + 

—  sin  Q  sin  a) 


is  given  by 

-  sin  Q  sin  cos  I 


—  sin  Q  cos  u)  cos  I 


cos  Q  sin  uu  cos  I 


+  cos  cos  0)  cos  I 
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Fig.  2.  Euler  angles  I,  fi ,  u; . 
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b23  = 

—  cos  n  sin  I 

b31  = 

sin  uu  sin  I 

b32 

cos  uu  sin  I 

b33  = 

cos  I 

(31) 

12  3 

(See  Ref.  3,  p.  328.  )  If  (y  ,  y  ,  y  )  is  a  point  in  the  orbit  of  the  body, 

3 

y  =  0  and 

1 

y  =  r  cos  i|f 

y^  =  r  sin  \|f  (3  2) 


By  (19)  -  (25)  and  (29)  we  have 


1  +  e  cos  ijf 


(33) 


(34) 

(35) 


III  -A. 


Parabolic  Motion 

In  the  case  of  parabolic  motion, 


1/2  -3/2 
M  p 


/d‘  -/■ 


(1  +  cos  i|r) 


e  =  1,  a  =  00  and  by  (34)  and  (35) 

=  i  (tan  I  +  r ta"3  l) 
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Let  t  be  the  time  of  perigee  passage  when  Jj 
P 

motion  n  and  mean  anomaly  M  at  time  t  by 


1/2  -  3/2 

n  =  U  P 


M  =  n(t  —  t  ) 
P 


=  0  .  We  define  the  mean 

(36) 

(37) 


Then  the  true  anomaly  i|/  at  time  t  is  determined  by  solving  the  equation 

1  *  ,3  i 


2M  =  tan  +  J  tan"  2 


as  follows  (see  Ref.  4,  p.  26).  Consider  the  identity 

W  ~^=  + H-r)3 


with 


\|r  ,  1 

tan  ^  -  X  —  — 


6M  =  X^  —  — 


Let 


X  =  —  tan  y 


X  =  —  tan  (3 


Then 


(3  =  cot  -1  (3M) 


Y  =  tan 


1  [(tan(3)l/3J 


-1 


^  =  2  tan  (2  cot  2y) 


0  *  e  *  ! 


—  T  ^  Y  5  ?" 
2  '  2 


—  IT  ^  l|f  ^  TT 


(38) 


(39) 
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Given  the  orbital  elements  (p,  e  =  1,  I,  0,  cju,  t  )  we  determine 

1  2  3  .1.2.3  P 

the  position  and  velocity  (x  ,  x  9  x  9  x  ,  x  9  k  )  at  time  t  by  first 

1  2 

determining  i|f  from  (37)  and  (39).  We  determine  r  by  (35)  and  (y  ,  y  ) 

3 

by  (32).  We  then  apply  (30)  with  y  =  0  to  determine  the  position 

coordinates.  The  velocity  is  given  by 


rJ  - 


-  E 


k=l 


b..y 

Jk 


1,2,3 


(40) 


where 


dr 


di|f 


=  gj-  cos  ♦  -  r  Sin  *  gf- 


dr  ... 

=  rr—  sm  Ur  +  r  cos 
dt 


d^ 

dt 


(41) 


with  dr/dt  and  dijf/dt  being  determined  by  (33)  and  (34).  dr/dt  <  0  if 

t  <  t  and  dr/dt  >  0  if  t  >  t 
P  P 

I1I-B.  Elliptic  Motion 

We  now  suppose  that  0  ^  e  <  1  and  a  >  0  .  Equations  (29)  and 
(3  0)  imply 


± 


where  the  mean  motion  n  is  defined  by 


n 


a 


-  3/2 


(42) 


(43) 


Let  t  be  the  time  of  perigee  passage  and  let  the  mean  anomaly  M  at 
P 

time  t  be  given  by 


M  -  n(t  —  t  ) 
P 


(44) 
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In  terms  of  the  mean  anomaly  M  at  time  t  we  have 

J  o  o 

M  =  M  +  n(t  -  t  ) 
o  o 

We  define  the  eccentric  anomaly  §  by 

r  =  a(I  -  e  cos 

with  ?  =  0  when  t  =  t  .  We  have 


P 

dr  =  a  e  sin  §  d  £ 


and  (42)  becomes 


M  =  n(t  -  t  )  -  ^  —  e  sin  5 

The  ±  sign  is  eliminated  by  changing  the  sign  of  £  if  necessary 
changing  (46).  See  Ref.  p,  22. 

Equations  (35  )  and  (45)  give 


.  cos  E  -  e 

COS  W  =  - - 2 - ; 

1  -  e  cos  c 


which  implies 


1  +  cos  tjj 


(1  -  e)  (1  +  cos  g ) 


1  -  e  cos  5 

1  ,  (1  +  e)  (1  -  cos  5) 

1  —  cos  ^ 

1  -  e  cos  ? 


These  equations  may  be  written 


2  cos 


2  ilf  1  -  e  0  2  § 

o  -  1 - 5"  •  2  cos 

2  1  -  e  cos  §  2 


9  •  2  f 

2  sin  £ 


1  +  e  0.25 

1 - =-  •  2  sm 

1  -  e  cos  5  2 


(45) 


(46) 


(47) 

without 


(48) 
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Dividing  the  second  equation  by  the  first  and  taking  square  roots  we  obtain 


tan  I  =Vf^i'  tan  I  (49) 

We  have  M  =  §  =  i|r  when  \|r  =  0  +  2krr  or  i|r  =  tt  +  2kn  ,  k  any 
integer.  Further,  modulo  2  tt  we  have 

0  <  ijf  <  tt  implies  0  <  ^  <  tt  ,  0  <  M  <  tt 

tt  <  i|r  <  2tt  implies  tt  <  §  <  2  tt  ,  tt  <  M  <  2  tt  (5  0) 


which  removes  any  quadrant  ambiguity  in  determining  ijf  from  ^  using 
(49).  The  body  completes  one  orbit  when  M  increases  by  2  tt  ,  so  that 


orbital  period  = 

Equations  (32),  (46)  and  (48)  imply 

yl  =  r  cos  i|j  =  a(cos  § 


2  tt 
n 


-  e) 


(51) 


(52) 


By  (46)  and  (52) 


2.2,  2 

r  s  in  ijr  =  r 


2  2 
r  cos 


/  2  2  2 
(1  —  2  e  cos  §  +  e  cos  §  —  cos  £ 


+  2  e  cos  £ 


-e2) 


2,i  2,  .  2 

a  (1  —  e  )  sm  £ 


so  that 


2 

y  =  r  sm 


e  sin  ^ 


because  of  the  relations  (50). 


(53) 
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Given  the  orbital  elements  (p,  0  ^  e  <  1  ,  I,  Cl,  u),  t  )  we 

1  2  3  •  1  •  2  •  3  P 

determine  the  position  and  velocity  (x,x,x,x,x,x)  at  time  t 
by  first  determining  M  from  (44)  and  reducing  it  to  be  between  0  and  2  n - 
Then  Kepler's  equation  (47)  is  solved  for  §  using  a  Newton-Raphson  itera¬ 
tion.  Namely,  if  f(§)  is  a  function  and  we  wish  to  find  §  such  that 
f(§)  =  M  starting  from  a  nearby  value  § Q  >  we  would  define  as  in  Fig.  3 

M  =  f(§  ) 

o  o 

M  -  M 

«o  * 

?,  =  +  A?0  (54) 

continuing  the  Newton-Raphson  iteration  to  define  M  ,  A§^  ,  »  etc., 

until  we  converge  to  §  .  Convergence  will  occur  unless  f7/(§)  =  0  •  In 

the  specific  case  of  Kepler's  equation  (47)  we  take  as  a  first  approximation 

5  =  M 

^  o 

If  e  =  0  ,  nothing  further  is  required.  Otherwise,  as  in  Ref.  5,  p.  84, 

we  let 


M  =  ?  —  e  sin  £ 

o  o  ^o 

M  -  M 
_ o__ 

”  1  -  e  cos  £ 

^o 

=  §0  +  A§q  (modulo  2tt)  (55) 

continuing  the  iteration  to  get  ,  A£^  t  §2  ’  etc*  *  until  we  obtain  a 

value  B  ,  such  that  I  A?,  I  <  €  ,  where  €  is  an  accuracy  constant 
^k+1  1  ^k  1  7 

depending  on  the  number  of  places  in  floating  point  computations  on  the 
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Fig.  3.  Newton- Raphs on  technique. 
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computer  being  used.  We  have 


t'  '  (§)  =  e  sin  § 


which  is  zero  only  at  §  =  0  and  §  =  tt  •  Thus  the  Newton-Raphson 

iteration  will  converge  everywhere  except  at  these  two  points.  But  the 
iteration  is  not  needed  at  these  two  points,  because  the  first  guess  §  =  M 

already  gives  the  correct  value  of  §  . 

1  2 

Having  obtained  £  we  evaluate  (y  ,  y  )  from  (52)  and  (53)  and 

12  3  3 

determine  (x  ,  x  ,  x  )  from  (30)  with  y  =  0  . 

Differentiating  (47)  with  respect  to  t  we  obtain 

“  <56> 


so  that  differentiating  (52)  and  (53)  we  obtain 


•  1 

y 


sin  l 


(57) 


2  na^  \fl  -  e 


T1 


cos  ^ 


(58) 


•  1  -2  *3 

The  velocity  (x  ,  x  ,  x  )  is  then  determined  by  (40). 


Ill  -  C .  Hyperbolic  Motion 

We  now  suppose  that  e  >  1  and  a  <  0  .  Equations  (27)  and  (33) 


imply 


n  |  a  |  ydt  =  J 


rdr 


y<r  +  I  a  I  )2  -  |  a  |  2  e 


(59) 
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where  the  mean  motion  n  is  defined  by 


1/2  ,  |-3/2 

n  =  4  a 


Let  t  be  the  time  of  perigee  passage  and  let  the  mean  anomaly 
P 

be  given  by 


M  =  n(t  —  t  ) 

P 

We  define  the  eccentric  anomaly  ^  by 

r  -  j  a  |  (e  cosh  §  —  1 ) 

with  F  =  0  when  t  =  t  .  We  have 

P 

dr  =  |  a  |  e  sinh  ^  d  § 


and  (59)  becomes 


M  =  n  (t  —  t  ) 
p 


e  sinh  5  5 


The  ±  sign  is  eliminated  by  changing  the  sign  of  §  if  necessary 
changing  (62).  See  Ref.  4,  p,  27. 

Equations  (27),  (35)  and  (67)  give 


cos  i|r 


e  -  cosh  g 

e  cosh  §  -  1 


which  implies 


1  4  cos 


(e  -  1 )  (1  4  cosh  l) 
e  cosh  §  -  1 


1 


COS  ijf 


(e  4  1 )  (cosh  g  -  1  ) 
e  cosh  §  -  1 


(60) 

M  at  time  t 


(61) 


(62) 


(63) 

without 


(64) 
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These  equations  may  be  written 


o  2  *  e  -  1  0  ,2? 

2  cos  ir  =  - r— - sr-  •  2  cosh  *- 

2  e  cosh  ^  -  1  2 


0  .  2  |  e  +  1 

2  sin  =  - 7— - r- 

2  e  cosh  ^  -  1 


•  2  sinh^ 


Dividing  the  second  equation  by  the  first  and  taking  square  roots  we  obtain 


tan 


je  +  1' 

'  Ve  -  1 


tanh  ^ 


(65) 


We  have  M  =  £  -  ^  when  \|r  =  0  .  Further 

0  <  ijr  <  tt  implies  0  <  ^  <  00  >  0  <  M  <  00 

—  n  <  f  <  0  implies  —  00  <  ^  <  0  ,  —  00  <  M  <  0  (66) 

which  removes  any  quadrant  ambiguity  in  determining  \(r  from  §  using  (65). 
Equations  (32),  (62)  and  (64)  implies 


y*  =  r  cos  i|r  =  |a|(e  —  cosh  5) 


(67) 


By  (62)  and  (67) 

2.2,  2  2  2,  II2/2  r2  o  2 

r  sm  ^  =  r  —  r  cos  ur  =  a  (e  cosh  ^  —  2e  cosh^  +  1  —  e 


+  2e  cosh  ^  —  cosh 

,,2  2  2 

=  a  ( e  -  1 )  sinh  £ 


so  that 


y  =  r  sm 


Ve2  -  1*  si 


inh  ? 


(68) 
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because  of  the  relations  (66). 

Given  the  orbital  elements  (p,  e  >  1 ,  I,  0,  r^,  t^)  we  determine 
the  position  and  velocity  (x  ,  x  ,  x  ,  x  ,  x  ,  x  )  at  time  t  by  first  deter¬ 
mining  M  from  (61).  Then  (63)  is  solved  for  §  using  a  Newton-Raphson 
iteration  (54).  Namely,  we  take  as  a  first  approximation 


=  M  if  |M|  £  1 
=  sinh  *(M/e)  if  |M|  >  1 

Then  we  let 

M  =  e  sinh  £  —  £ 

o  ~  o  - 

M  -  M 

A?  3  - TT - 2T“ 

~ o  e  cosh  £  -  I 

5  o 

§1  =  5o  +  «o  (69) 


continuing  the  iteration  to  get  ,  A^j 
Z,  ,  ,  such  that 


2’ 


etc.  ,  until  we  obtain  a  value 


«ki <  £  if  i5k+ii s  1 


|A?kl 

I -k+1 I 


if 


>  1 


where  e  is  an  accuracy  constant  depending  on  the  number  of  places  in 
floating  point  computation  on  the  computer  being  used.  We  have 


f'  7  (?)  =  e  sinh  § 
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which  is  zero  only  at  5  =  0  •  Thus  the  Newton-Raphson  iteration  will  con¬ 

verge  everywhere  except  at  this  point.  But  the  iteration  is  not  needed  at 
this  point,  because  the  first  guess  =  M  already  gives  the  correct  value 

of  e;  . 

1  2 

Having  obtained  5,  we  evaluate  (y  ,  y  )  from  (67)  and  (68)  and 

12  3  3 

determine  (x  ,  x  ,  x  )  from  (3  0)  with  y  =  0  . 

Differentiating  (63)  with  respect  to  t  we  obtain 


d£  _  n | a | 
dt  r 


so  that  differentiating  (67)  and  (68)  we  obtain 

,  .2 


.  1 

y 


n  a 


sinh  § 


.2 

y 


cosh  ^ 


The  velocity 


is  then  determined  by  (40). 


(70) 


(71) 


(72) 


III-D.  Straight  Line  Motion 

In  the  case  that  G  =  "5  (see  equation  (13)),  the  vectors  r  and  v 
are  always  parallel.  Thus  we  have 

-  dr  r 
V  ~  dt  r 


and  (12)  implies 


v  = 


dr 

dt 


±V2<  r  +  C)' 


(73) 
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so  that  if  the  initial  time  is  t 

o 


t  -  t 

o 


=  /  dt  =  /  ;  dr 

+  C)' 


3C2r2  —  4C|jr  + 

3 

15C 


+  Cr‘  +  c. 


(74) 


if  the  velocity  does  not  change  sign  between  tQ  and  t  •  Let  the  initial 
position  vector  r  have  components 


x  =  r  cos  d  cos  9 
o  o 


2 

x  =  r  cos  d  sin  0 
o  o 


3 

x  =  r  sin  rt 
o  o 


(75) 


Let  the  initial  velocity  vector  have  magnitude  v 


Then  by  (73) 


C 


2  Vo 


u 

r 

o 


(76) 


and  by  (73) 

3C^r  ^  —  4Cfir  +  8(j^ 

c,  =  +  - 2 - 5 - 2 -  yTTUTTcP  (77) 

1 5CT  y 


Let  us  take  as  orbital  elements  in  the  straight  line  case 

(r  ,  e  =  oo  Q  M  v  ,  t  ).  Then  using  the  above  formulas  we  can 

o  o  o 

1  2  3  . 1  . 2  •  3 

determined  (x,x,x,x,x,x)  at  time  t  from  the  orbital  elements. 
Conversely  it  is  trivial  to  determine  the  orbital  elements  from  the  position 
and  velocity  at  time  t 
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The  above  formulas  should  be  put  in  better  form  for  computations, 
but  since  we  will  really  never  be  concerned  with  the  straight  line  case,  we 
shall  not  bother  to  do  so. 


IV.  DETERMINATION  OF  CONIC  SECTION  ORBITAL 

ELEMENTS  GIVEN  THE  POSITION  AND  VELOCITY 


Suppose  we  are  given  the  position  and  velocity  (r^,  v^)  = 

Z  3*1  2  •  3 

x  ,  x  ,  x  ,  x  ,  x  )  at  time  t  and  wish  to  determine  the 
o  o  o 


(X  ,  -Tk.  •  ^  f 

o  o  o  o  9 


conic  section  elements  (p,  e,  I,  Q,  (u,  t  )  as  defined  in  (2)  and  (3). 

We  first  calculate  the  angular  momentum  vector  as  given  in  (13), 


G  =  r  xv 
o  o 


(78) 


If  G  =  |fi|  =  0  ,  we  have  the  straight  line  case  for  which  the  determina¬ 

tion  of  the  elements  follows  easily  from  the  discussion  in  Section  III-D. 

If  G  /  0,  we  determine  the  semi-latus  rectum  p  from  (22). 

By  (29)  we  have 


1 


a 


(79) 


and  by  (27)  we  have 


e  = 


(80) 


The  inclination  I  and  ascending  node  Q  are  determined  by  (15) 
and  (16).  If  I  =  0°  or  180°,  fl  is  indeterminant  so  we  can  take  fi  =  0 

by  convention  in  this  case. 
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Solving  (30)  and  (40)  for  b  f  we  have 

jk 


j2 


1  /  j.2  . i  2' 

— r."  2 — n— r  (x  y  ~  x  y 

y0  yo  ~  yG  yo  v 


rr"2  n  z 

y  y  -  y  y 

o  o  ^o  o 


j  =  1,2,3  (81) 


•j  1  j*1 

x  y  —  x  y 


By  (31),  we  have  for  j  =  3  , 


sin  cu  = 


cos  uu 


31 


sin  I 
b„ 


32 


sin  I 


(82) 


If  I  =  0°  or  180°  so  that  by  convention  fl  =  0  ,  we  can  use 


sin  uu 


12 


=  ±  b 
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cosw  =  bn  =  ±  b22 


(83) 


where  the  +  sign  of  the  ±  sign  is  to  be  used  if  I 
to  be  used  if  I  =  180°  . 


0  and  the  —  sign  is 


In  order  to  determine  the  time  of  perigee  crossing  t  and  the  com- 

1  2  .1  .2 

ponents  of  position  and  velocity  in  the  orbital  plane  y  ,  y  ,  y  ,  y  for 
determining  the  argument  of  perigee  ud  from  (81),  (82)  and  (83),  we  must 
consider  three  separate  cases. 


IV -  A .  Parabolic  Motion 

If  e  =  1  or  if  !  e  -  1  j  <  e  ,  where  e  is  an  accuracy  constant 
depending  on  the  number  of  places  in  floating  point  computations  on  the 
computer  being  used,  then  we  are  in  the  parabolic  case  and  would  set 
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e 


1  identically.  Equation  (35)  is 


r 


_ E _ 

1  +  COS  l|f 


(84) 


and  (33)  with  a  =  00  is 

raf  =  v2w  -  w 

Differentiating  (84)  and  using  (34)  we  have 


dr 

dt 


sin  tjr 


We  thus  have  at  the  initial  time  t 

o 


cos  i j 
T  o 


sin  ilf 

o 


V2Pro  -  P2  ' 

r 

o 


(85) 


from  which  li;  is  easily  determined  with  —  tt  <  ilf  <  tt 
o  J  o 

the  time  t  of  perigee  passage  is  then 


t 

P 


t 

o 


+ 


By  (37)  and  (38) 


(86) 


1  Z  •  1  •  2 

where  n  is  given  by  (36).  The  coordinates  y^  ,  yQ  ,  y^  ,  yQ  are  given 
by  (32)  and  (41)  for  use  in  determining  cjo  ♦ 


IV “B.  Elliptic  Motion 

We  now  suppose  that  0  ^  e  <  1  and  a  >  0  .  If  e  <  e  ,  where 
e  is  an  accuracy  constant,  we  set  e  =  0  identically  for  the  circular  case. 

If  e  =  0  ,  we  would  take  t  =  tQ  and  for  use  in  determining  uu 

assume  that 
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1 


=  0 


.  1 

y„ 


=  o 


y  =  v 


If  e  ^  0,  we  differentiate  (46)  with  respect  to  time  to  obtain 


dr  .  _ 

ar  =  sin  5  at 


d  g  na  e  sin  g 
r 


by  (56),  Since 


dr 

dt 


r  .  v 


this  gives  at  the  initial  time 

sin  §  = 

Further,  (46)  can  be  put  in  the  form 


r  •  v 
o  o 

Z 

na  e 


COS  5  =  —  (1 - — ) 

e  a  ' 


Then  the  time  t  of  perigee  passage  is  given  by  (47)  in  the  form 
P 


t  =  t  —  —  (£  —  e  sin  £) 

p  o  n 


1  2  1*2 

and  yQ  ,  y^  ,  yo  ,  y^  are  given  by  (52),  (53),  (57),  (58)  for  use  in 
mining  u)  . 


(87) 

(88) 

(89) 

(90) 

(91) 

(92) 

deter  - 
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IV -C.  Hyperbolic  Motion 


We  now  suppose  that  e  >  1  and  a  <  0  .  Differentiating  (62)  with 
respect  to  time  we  obtain 

2 

dr  (|  .  1  r  n  |  a  |  e  sinh  § 

dF  =  la|  e  sinh  5  af  =  -  (93) 

by  (70)  .  Using  (89),  we  have  at  the  initial  time 

r  •  v 

sinh  §  = 

■  n|a|2e 

Further,  (62)  can  be  p\it  in  the  form 

,  r 

cosh  §  =  —  (1  +  “2-) 

*  e  a 

Then  the  time  t  of  perigee  passage  is  given  by  (63)  in  the  form 

t  =  t  —  —  (e  sinh  §  —  £)  (96) 

p  o  n  v  °  •  '  ' 

1  2  •  1  •  2 

and  y^  ,  y^  ,  y^  ,  y^  are  given  by  (67),  (68),  (71),  (72)  for  use  in  deter¬ 
mining  (jo  . 

V.  LAMBERT  PROBLEM  FOR  THE  TIME  OF  TRAVEL 

BETWEEN  TWO  POINTS  IN  A  CONIC  SECTION  ORBIT 

Lambert’s  problem  is  to  determine  the  elements  (p,  e,  I,  Q,  oo,  t  ) 

— ♦  P 

of  the  conic  section  orbit  which  has  position  vectors  r  ^  at  time  t^  and 
r?  at  time  t^  with  t^  <  t^  •  We  shall  consider  a  variation  on  the  prob¬ 
lem  in  that  we  shall  take 

h  ’  V1  =  I  Vj  I  ,  tj  ,  r2  (97) 

as  given  and  determine  the  elements  (p,  e,  I,  fi,  uo,  t  )  and  the  time  of 

P 

flight  ( 1 2  —  tj)  from  to  r^  . 


(94) 


(95) 
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If  r*2  and  are  parallel  (r^  =  kr^  with  k  ^  0),  then  we  have 
the  case  of  straight  line  motion  with  elements  easily  determined  by  the 
formulas  in  section  III-D.  If  r^  and  r^  are  anti-parallel  (r  =  kr^ 
with  k  <  0),  then  the  orbital  plane  is  indeterminant  and  there  are  infinitely 
many  orbits  (all  with  the  same  elements  p,  e,  t^,  t^  —  tj )  satisfying  the 
conditions  (97), 

Now  let  us  suppose  r^  x  r^  ^  "3*  The  angle  9  from  r  to  rn 
rotating  in  the  direction  of  motion  is  then  determined  by 


sin  0  =  ± 


cos  0  = 


rl  X  r2 


(98) 


where  the  +  sign  is  used  if  0  <  0  <  2tt  and  the  —  sign  is  used  if 
tt  <  0  <  2tt  .  Since  we  are  concerned  with  rapid  intercept  orbits,  we  shall 
assume  0  <  0  <  n  in  the  program,  although  for  the  sake  of  completeness 
we  shall  allow  both  cases  in  the  formulas  that  follow.  The  distance  c 
between  the  points  r^  and  r^  is  given  by 


2  2 

*  rl  +  r2  -  2rl 


(99) 


The  unit  vector  normal  to  the  orbital  plane  is 


A 

G  = 


(100) 


where  the  +  sign  is  used  if  0  <  0  <  tt  and  the  —  sign  is  used  if 

tt  <  0  <  2tt  .  The  inclination  I  and  ascending  node  Q  are  determined 

A 

from  G  by  formulas  analogeous  to  (15)  and  (16). 


3  1 


By  (29),  the  semi-major  axis  a  is 


1_ 

a 


and  the  unknown  magnitude  of  velocity 

Z 

v2  =  M 


is 


(101) 


(102) 


If  the  right  side  of  (102)  is  negative  or  zero,  a  conic  section  orbit  satisfying 
(97)  does  not  exist.  By  (29)  and  (33)  we  have 


+  “^2  =  Vj2  ,  i  =  1,2  (103) 

r . 

1 

Now  consider  unit  vectors  ~e  ^  in  the  orbital  plane  with  ~e  ^ 

pointing  towards  perigee  and  with  normal  to  e ^  pointing  in  the  direc¬ 
tion  of  motion  at  perigee.  We  have 

r.  =  y.1  Sj  +  y?  e2  ,  i  =  1,2  (104) 


x^  )  are  the  components  of  r^  in  the  usual  coordinate  system, 


Xij  =  t  b..y.k  ,  i  =  1,2  ;  j  =  1,2,3  (105) 

k- 1  J 


Solving  for  b  ^  ,  we  obtain 


h  _  _ 1 _ 

jl  "  12  12 

Yi  y2  -  y2  y i 


v  _  i 

j2  12  12 

yi  y2  -  yz  yi 


/  j  2  j  2^ 

X1  y2  "  x2  yl 


j.  1  -  .  L  lN 


•x2  yi  -  vy2 


j  =  1,2,3 


(106) 


If  (x.1,  x.2, 

(30)  gives 
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The  argument  of  perigee  (*)  could  then  be  determined  by  (82)  and  (83)  if 
y^  ,  y?  were  known. 

To  recapitulate,  we  are  given  values  of  the  parameters  (97)  from 
which  we  determined  I  ,  Q,  a  and  v^  .  In  order  to  determine  (p,  e,  uo, 
tp,  t^)  we  shall  have  to  consider  separately  the  cases  of  parabolic  motion 
(a  =  oo  )f  elliptic  motion  (a  >  0)  and  hyperbolic  motion  (a  <  0)  making 
use  of  relations  (98)  through  (106). 

V  -A.  Parabolic  Motion 

If  a  =  oo  so  that  e  =  1,  we  have  by  (35) 


ri  1  +  cos  . 


p  £  i 
2SeC  ~ 


i  =  1,2 


(107) 


Then  by  (34)  and  (107) 


r^  =  1  sin  \J l  ,  i  =  1 , 


(108) 


so  that  by  (107)  and  (108) 


cos 


sm 


=  JL  -  1 


r . 

l 


= 

U  1 


i  =  1,2 


(109) 


From  (37)  and  (38)  we  obtain 


1  /  *i  1  3  +i\ 

n(ti  “  tp)  =  2  \tanT  +  3  tan  ~) 


,  i  =  1,2  (110) 


with 


0  =  to-  ~ 


(111) 
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By  (32),  equations  (104)  take  the  form 


+  — *  f  ~ * 

=  r.  cos  i t.  e,  +  r.  sm  i e0  ,  1  =  1,2 

l  l  Yi  1  i  Yi  2 


so  that  by  (99) 


c  =  Tj  +  r2  “  2rir2  cos (^2  ~  ^  1  ^ 


2  2r2  v  1\ 

=  (rj  +  n2)  —  2v^Q.^  cos  ( — - I 


which  implies 


■vFl‘2  \  2 


2  A/r ,  r  *  cos 


*2  -  *i\ 


=  +  r2  +  c)  (FJ  +  r2  -  c)1  (112) 


where  the  +  sign  is  used  if  0  =  iff ^  ^ j  <  TT  and  the  —  sign  is  used  if 

0  =  $2  ““  >  tt  .  Substituting  (107)  into  the  left  side  of  (112)  we  obtain 


tj  ^ 2  V( ri  +  r2  +  C)  (rl  +  r2  ‘  c) 

1  +  tan  — r—  tan  =  ±  1 - 

2  2  p 


(113) 


It  also  follows  from  (107)  that 


,  E  h  x  f  2  1  X  «.  2  21 

rl  +  r2  =  2  2  +  tan  T  +  tan  T 


The  last  two  equations  give 


<rl  +  r2  +  c)  +  (rj  +  r2 


—  c)  T  +  r2  +  c)  (rx  +  -  c) 1 


=  (tan  — =r  —  tan 
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which  implies 


+  r2  +  c'  TVrl  +  r2  ~  C'  *2  h 

-  =  tan  —  tan  -g- 

(114) 

Subtracting  the  equation  of  (1  1 0)  with  i  =  1  from  the  one  with  i  =  2  , 

we  obtain 


6n(t2  ~  tj)  = 


/  tjf  2  1N 

ltanT  “  tan~ 


/  *1  M 

3  11  +  tan  —  tan  —I 


+  tan 


tan 


(115) 


Equations  (36),  (113),  (114)  and  (115)  yield 


6|j  ^2(t2  ~~  ti )  =  (ri  +  r2  +  c)3//  T  (r^  +  r2  —  c)3/2  (116) 


where  the  —  sign  is  used  if  0  <  0  <  tt  and  the  +  sign  is  used  if 
rr  <  0  <  2jt  .  This  result,  determining  t^  in  terms  of  t^  ,  is  called 
Euler's  equation  for  parabolic  motion  and  is  given  in  Ref.  3,  pp.  157-158. 
By  (1  09)  we  have 


so  that  (114)  becomes 


sin  tlf . 

_  i 

1  I  cos  f. 

1 


r .  r . 
i  i 


(117) 


r2r2 


rlrl 


(118) 
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Equations  (111),  (113),  (117)  and  (118)  yield 


*2  *1 
1  +  tan  —  tan  — 


rl  +  r2  +  c 


-4 


rl  +  r2  -  C 


yf( rj  +  r2  +  c)  (r^  +  -  c)' 


(H9) 


where  the  +  sign  of  the  ±  is  used  if  0  <  9  <  tt  and  the  _  sign  is  used  if 

n  <  0  <  2tt  .  All  the  quantities  in  (119)  except  p  are  known,  so  that  (119) 
determines  p  .  Then  r^  ,  r^  are  determined  by  (103)  in  the  form 


*i  -  «*>i 

r.  i  r. 


=  1.2 


(120) 


where  the  +  sign  is  used  if  t.  <  t  and  the  —  sign  is  used  if  t.  >  t 

1  P  1  P 

However,  whether  T  is  before  or  after  the  time  of  perigee  t  is  not  yet 


known,  so  r\  is  determined  only  up  to  sign.  The  tangent  of  is  given 

by  (117)  up  to  sign  and  the  value  of  (t^  —  t  )  is  given  by  (110)  up  to  sign. 
We  know  that 


(±)2 


^2  1  3  ^2 

"  Wi 

h 

i 

3  h 

tan  +  —  tan  -g- 

tan  — 

+  3 

tan  ^ 

2n(t2  -  tx) 

(121) 


where  the  right  side  is  known  from  (116).  In  general  there  is  only  one 
combination  of  signs  which  will  make  (121)  valid,  so  we  may  regard  the 

♦  i 

sign  ambiguity  for  r^  and  tan  as  resolved.  Then  t  is  determined  by 
(110),  cos  i|k  and  sin  by  (109),  yj  by  (32)  and  uu  by  (106),  (82)  and  (83). 
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V  -B.  Elliptic  Motion 

If  a  >  0  so  that  0  ^  e  <  1  we  have  by  (46) 


=  a(l  —  e  cos  §.)  ,  i  =  1 ,  2 


Then  by  (56)  and  (122) 


r .  = 

1 


yJ]IP  e  sin 


r . 
i 


,  i  =  1,2 


Equations  (122)  and  (123)  can  be  written 


r . 


e  cos  F .  =  1  —  — — 

a 


e  s  in  §  • 


r .  r . 
i  i 


1 

Squaring  these  equations  and  adding  we  obtain 

2 


2  2.2 

r.  \  r.  r. 

II - L]  +  — — 

a  /  |ia 


Kepler's  equation  (47)  implies 


i  =  1,2 


i  =  1,2 


(122) 


(123) 


(124) 


(125) 


n(t.  —  t  )  =  F.  —  e  sin  F .  ,  i  =  1,2 

l  p  5i  -i 


By  (52)  and  (53)  equations  (104)  take  the  form 
r.  =  a(cos  F.  —  e)  c,  + 

l  l  1 


(126) 


- 2? 

1—  e  sin  F.  e.  ,  i  =  2 

2 
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so  that  by  (99) 


2  2  2  2  2.  2 
c  =  a  (cos  §2  “  cos  §j)  +  a  (1  —  e  )  (sin  “  sin 


=  4a 


2  jj _ 2 _ 2  (?l  +  ^ 


e  cos 


2  (?2  ?l) 


sm 


It  follows  from  (122)  that 


rl  +  r2  *  2a 


1  —  e  cos 


(Si  +  ?*> 


cos 


we  obtain 


n(t2  -  tj)  =  (§2  -  fj)  - 


2  e  cos 


2 

^2 

-  «!>' 

2 

the 

one  with  i  =  2 

+  ?2> 

.  (-2 

(127) 


(128) 


(129) 


Continuing  to  follow  Ref.  7,  we  define  the  angles  a  and  (3  by 

(?1  +  S2) 


cos 


(a  +  P) 


=  e  cos 


0  ^  a  +  (3  <  2tt 


a  —  (3  =  —  §i  —  2mrr  0  ^  a  —  p  <  2tt 


(130) 


where  m  is  the  number  of  complete  circuits  made  between  times  t^  and  t2 
Inequalities  (13  0)  imply 


0  ^  a  <  2tt  ,  —  tt  ^  P  <  tt 

Equations  (127),  (128)  and  (129)  become 


(131) 


c  .  a  +  P  .  a  -  p 

tz —  —  Sin  1 ^  sm  1 - - — ^-L 

2a  2  2 


(132) 


38 


1 


cos 


(133) 


rl  +  r2 
2a 


—  cos 


(Q.  +-.P-)  roR  (a  -  P) 


)  =  2mrr  +  a  —  (3  —  2  cos  — — ^  ^  sin  — — - — — 


(134) 


There  is  no  ambiguous  sign  in  (132)  because  by  (130) 


0  s  i-i-S-  <  n  ,  sin  (2-1-21  *  0 


0  s  °  2  P  <  n 


,  sin  !“  .  T  fij  *  o 
2 


(135) 


Equations  (132)  and  (133)  imply 


sm  2 


,'Vrr 


+  r2  +  c' 


2V&1 


(136) 


.  P 

sm2  = 


V 


ri  + 


—  C 


2^a1 


(137) 


so  that  a  and  (3  are  determined  in  terms  of  the  known  quantities  r^,  r^, 
c  and  a,  except  for  ambiguity  in  sign  in  (137).  There  is  no  ambiguity 
in  (136)  because  of  (131).  If  either  the  right  side  in  (136)  or  that  in  (137) 
are  greater  than  1  in  absolute  value,  then  an  elliptic  orbit  passing  through 
the  two  given  points  with  the  given  magnitude  of  velocity  at  the  first  one 
does  not  exist.  Equation  (134)  can  be  written  in  the  form 


n(t^  —  tj )  =  2mn  +  (a  —  (3)  —  (sin  a  —  sin  (3)  (138) 

which  is  called  Lambert's  theorem  in  Ref.  4,  p.  51.  We  have  thus 
determined  the  time  t^  except  for  the  number  m  of  revolutions  and  the 
ambiguity  of  sign  in  (137).  For  the  rapid  intercept  problem  of  interest  to 
us  we  can  assume  m  =  0  . 
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By  (111)  and  (49)  we  have 


tan  2"  =  tan 


'*2  ~  *1 


f  2  tj 

tan  -y  -  tan  — 

1  +  tan  —  tan  — 


i — V  r  ^  2  ?r 

V 1  -  e  tan  -y  -  tan  -y 


(1  -  e)  +  (1  +  e)  tan  —  tan 


h_ 

2 


yJiT 


T1 


-2  •  h  ?2 

sm  -y  cos  —  -  sm  —  cos  — 


^2  ^1  *2  .  ^1 
(1  -  e)  cos  -y  cos  -y  +  (1  +  e)  sin  -y  sin  — 


'1  - 


2  1 

e  sin 


'§9  "  § 


COS 


?2  "  -11 


+  §• 


-  e  cos 


By  (130)  this  becomes 


tanf 


xr  v  . 

-  e  sm 


cos 


m 

-  “■  (H1) 

in(2-J) 


1  -  e  sin 


,  .  a  .  (3 

2  sin  y  sin  y 


(139) 


Thus  by  (135),  (136)  and  (139)  we  can  resolve  the  sign  ambiguity  in  (137) 
as  follows: 

(i)  0  <  0  <  TT  -  sin  y  >  0 
(ii)  TT  <  0  <  2tt  -*  sin  y  <  0 
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We  now  calculate  the  remaining  elliptic  orbital  elements.  The 
eccentricity  e  is  determined  by  solving  (139),  obtaining 


2 

e 


20.2a  .  2(3 

tan  ^  sm  ^  sm  2 


sin 


m 


(140) 


The  semi-latus  rectum  p  is  given  by  (27),  The  quantities  r^  and  r^  are 
determined  by  solving  (125),  obtaining 


r.  =  (±)i^^Ve2  “  (!  -  “)  '  i  =  1,2  (141) 

i 

where  the  ambiguous  sign  (±)^  for  r.  is  to  be  resolved.  The  values  of 
sin  ^  (with  ambiguous  sign  (±)^)  and  cos  £ ^  are  given  by  (124)  if  e  >  0  . 
Subtracting  the  equation  of  (126)  with  i  =  1  from  the  one  with  i  =  2  and 

comparing  with  (130)  and  (138)  we  have 

(±)2|sin52|  -  (±)j  |  sin  §  j  |  =  Sm  ae~  Sm  ^  (142) 


In  general  there  is  only  one  combination  of  signs  which  will  make  (142) 

valid,  so  we  may  regard  the  sign  ambiguity  as  removed. 

Knowing  sin  and  cos  5^  (i  =  1 ,2),  we  can  determine  and 

^2  subject  to  the  restriction  provided  by  the  second  equation  of  (130)  with 

m  =  0  in  the  rapid  intercept  case.  Then  t  is  determined  by  either 

1  ^2 

equation  in  (126).  The  quantities  y^  and  y^  (i  =  1,2)  are  determined 

by  (52)  and  (53),  and  the  argument  of  perigee  au  is  determined  by  (106),  (82) 
and  (83). 

If  e  =  0  or  e  <  €  ,  where  e  is  an  accuracy  constant  depending 

on  the  number  of  places  in  floating  point  computations  on  the  computer  being 
used,  we  are  in  the  circular  case  with  r ^  and  would  set  e  =  0 

identically.  Since  there  is  no  perigee  point,  t  is  arbitrary,  so  we  can 
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W e  have 


take  t  =  t,  . 
p  1 


ri 

0 

1*1  cos  0 
Tj  sin  0 


(143) 


and  the  argument  of  perigee  co  is  determined  by  (106),  (82)  and  (83). 

V -C.  Hyperbolic  Motion 

If  a  <  0  so  that  e  >  1  we  have  by  (62) 


r.  =  |a  I  (e  cosh  —  1 )  ,  i  =  1 , 2 


(144) 


Then  by  (70)  and  (144) 


r . 

l 


^\Jl  |  a  I’  e  sinh 


r . 

l 


,  i  =  1,2 


(145) 


Equations  (144)  and  (145)  can  be  written 


e  cosh  E.  =  1  +  t~ r 

1  FT 


e  s  inh  £  • 


r .  r . 

l  l 


1  V7R 


i  =  1,2 


(146) 
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Squaring  these  equations  and  subtracting  the  second  from  the  first  we  obtain 


2 

e 


2.  2 

r .  r . 

l  l 

dial 


1,2  (147) 


Kepler's  equation  (63)  implies 


n(t.  —  t  )  =  e  sinh  ,  i  =  1,2  (148) 

By  (67)  and  (68)  equations  (104)  take  the  form 


cos 


+ 


1  sinh  £  e 


so  that  by  (99) 


1,2 

(149) 


2 

c 


|  a  [  ^(e^  —  l)  (sinh  ^  ““  sinh  §  j) 


(150) 


It  follows  from  (144)  that 


(151) 

Subtracting  the  equation  of  (148)  with  i  =  1  from  the  one  with  i  =  2  we 

obtain 


e  cosh 


?!  +  ?. 


cosh 


*2  -  ? 


-  1 


n(t2  “  h) 


2e  cosh 


sinh 


(152) 
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We  define  quantities  a  and  (3  by 


cosh 


(Hr6)  =  e  cosh  f 1 2  ~2) 


a  -  (3  =  -  § 


a  +  (3  >  0 

a  —  (3  >0 


(153) 


Inequalities  (153)  imply 


a  >  0 


(154) 


Equations  (150),  (151)  and  (152)  become 


-rpri 


=  s 


inli  sinh  (^1) 


(155) 


rl  +  r2 


TJH 


=  C 


°shH^)cosh(2x1)  - 


(156) 


n(t2  —  tj)  =  2  cosh  (—t  ^)  sinh  (a  ~  P)  —  (a  —  (3)  (157) 

There  is  no  ambiguous  sign  in  (155)  because  by  the  inequalities  in  (153), 
the  hyperbolic  sines  are  always  positive.  Equations  (155)  and  (156)  imply 


sinh  ^  - 


a  +yr,  +  r,  +  c' 


1  l2 


2/T 


(158) 


p  Vl  +  r2  -  c 

smh  p  =  - 


(159) 


44 


so  that  a  and  (3  are  determined  in  terms  of  the  known  quantities  ,  r^, 
c  and  a,  except  for  the  ambiguity  in  sign  in  (159)*  There  is  no  ambiguity 
in  (158)  because  of  (154).  Equation  (157)  can  then  be  written  in  the  form 


n(t^  ~~  tq)  =  (sinh  a  —  sinh  (3)  —  (a  —  (3) 


(160) 


We  have  thus  determined  the  time  t^  except  for  the  ambiguity  of  sign  in 

(159). 

By  (111)  and  (65)  we  have 


t  2  y  1 

6  .  *2-+l  tan  —  -  tan  — 

tan  2  =  tan  — j -  =  - j- - 

1  +  tan  —  tan  — 


1 1  •  ^ 


Je  -  1  Itanh  -  tanh 


a 


s  1 

(e  -  1 )  +  (e  +  1 )  tanh  -75-  tanh  — 


v, 


2 — T1 

e  -  1 


^2  ^1  bl  ^2 

sinh  cosh  — —  sinh  -g-  cosh 


^2  ^1  ^2  .  ^1 
(e  -  1 )  cosh  —  cosh  -j-  +  (e  +  1 )  sinh  -g-  sinh 

JJ7?  sinh(^-l±) 
ecosh(Eilli).cosh(^ 


By  (153)  this  becomes 


Ve2  -  l'  sinh (2-^) 

h(^-2)-  cosh(i^l) 


0  ¥  . v  2 

tan  j-  =  - 

cos 


inh(2-rfi) 


2  sinh  sinh 


(161) 
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Thus  by  the  inequalities  in  (153)  and  (154)  we  can  resolve  the  sign  ambiguity 
in  (159)  as  follows: 


(i)  0  <  0  <  TT  -*  sinh  —  >  0 
(ii)  tt  <  0  <  2tt  -*  sinh  —  <  0 


We  now  calculate  the  remaining  hyperbolic  orbital  elements.  The 
eccentricity  e  is  determined  by  solving  (161),  obtaining 


2 

e 


a  ,  2  0  ..2a  .  ,26 

4  tan  2"  sin^  sinh  — 


sinh' 


(162) 


The  semi-latus  rectum  p  is  given  by  (27).  The  quantities  r^  and  r^  are 
determined  by  solving  (147),  obtaining 


where  the  ambiguous  sign  (±)^  for  r^  is  to  be  resolved.  The  values  of 
sinh  (with  ambiguous  sign  (±)^)  and  cosh  are  given  by  (146).  Sub¬ 
tracting  the  equation  of  (148)  with  i  =  1  from  the  one  with  =  2  and 

comparing  with  (153)  and  (160)  we  have 


(±)2|sinh§2|  -  (±)j  |  sinh  §1  |  = 


sinh  a  -  sinh  ft 
e 


(164) 


In  general  there  is  only  one  combination  of  signs  which  will  make  (164) 

valid,  so  we  may  regard  the  sign  ambiguity  as  removed. 

Knowing  sinh  we  can  determine  §.  .  Then  t  is  determined  by 
1  l  1  2  p 

either  equation  in  (148).  The  quantities  y^  and  y^  (i  =  1,2)  are 

determined  by  (67)  and  (68)  and  the  argument  of  perigee  co  is  determined  by 
(106),  (82)  and  (83). 
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VI. 


LOGICAL  FLOW  OF  THE  COMPUTER  PROGRAM  TO  CALCULATE 
TIME  TO  INTERCEPT  VERSUS  INITIAL  VELOCITY 


1.  Initialize  constants. 

2.  Read  input  data  described  in  Section  I;  end  program  if  no 
more  data. 

3.  Calculate  position  and  velocity  of  launching  site  at  launch  time 
(earth  launch  uses  the  formulas  in  Section  II  and  parking  orbit  launch  uses 
the  formulas  in  Section  III-B.  ). 

4.  For  magnitude  of  velocity  v^  from  initial  magnitude  VLNCHO 
to  final  magnitude  VLNCHl  at  increments  DVLNCH,  calculate  and  store  in 
arrays  the  following: 

a.  Time  to  intercept  the  target  satellite  with  given  magnitude 
of  launch  velocity  v^  in  the  intercept  orbit  using  the  formulas  in  Section  V; 
calculated  iteratively  because  the  target  satellite  is  allowed  to  move  in  its 
orbit  between  launch  epoch  and  intercept  epoch. 

b.  Velocity  imparted  at  launch,  which  is  equal  to  the  vector 
velocity  (of  magnitude  v^ )  in  the  intercept  orbit  at  launch  calculated  using 
the  formulas  in  Section  III  minus  the  vector  velocity  of  the  launching  site. 

c.  Position  and  velocity  in  the  intercept  orbit  at  possible 
midcourse  correction  epochs  calculated  using  the  formulas  in  Section  III. 

5.  Plot  time  to  intercept  versus  magnitude  of  velocity  imparted  at 
launch.  A  point  is  deleted  from  the  plot  if  the  time  of  perigee  in  the  inter¬ 
cept  orbit  is  between  the  times  of  launch  and  intercept  and  if  the  perigee 
distance  is  less  than  the  radius  of  the  earth  plus,  say,  50  kilometers. 

6.  If  midcourse  corrections  are  to  be  considered,  change  the  target 
satellite  orbital  elements  by  the  input  orbital  element  errors. 

7.  For  each  original  intercept  orbit  calculated  in  (4.  )  and  each 
midcourse  correction  epoch  with  midcourse  launching  position  and  velocity 
calculated  in  (4c.  ),  do  the  following: 


47 


a.  For  magnitude  of  velocity  from  an  initial  to  a  final 

magnitude  at  given  increments  calculate  and  store  in  arrays  the  following: 

i.  Time  to  intercept  the  new  target  satellite  with 
given  magnitude  of  midcourse  launch  velocity  v^  in  the  intercept  orbit 
using  the  formulas  in  Section  V;  calculated  iteratively  because  the  target 
satellite  is  allowed  to  move  in  its  new  orbit  between  midcourse  correction 
epoch  and  intercept  epoch. 

ii.  Velocity  imparted  at  midcourse  correction  which 
is  equal  to  the  vector  velocity  (of  magnitude  v^ )  in  the  new  intercept  orbit 
at  the  midcourse  epoch  calculated  using  the  formulas  in  Section  V  minus 
the  vector  velocity  in  the  original  intercept  orbit  at  the  midcourse  epoch. 

b.  Plot  time  from  midcourse  correction  to  intercept  versus 
magnitude  of  velocity  imparted  at  midcourse  correction.  Points  are  deleted 
which,  by  the  criterion  described  in  (5.),  represent  intercept  orbits  inter¬ 
secting  the  earth’s  atmosphere. 

8.  Go  to  (2.  ). 
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